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Discussion of “The Potential and Limitations 
of Direct and Large-Eddy Simulations” 
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NASA Langley Research Center 
Hampton, VA 23665 

T. A. Zang 
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NASA Langley Research Center 
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Abstract 

The full text of the discussion paper presented at the Whither Turbulence Work- 
shop (Cornell University, March 22-24, 1989) on the potential and limitations of direct 
and large-eddy simulations is provided. Particular emphasis is placed on discussing 
the role of numerics and mathematical theory in direct simulations of both compress- 
ible and incompressible flows. A variety of unresolved issues with large-eddy simu- 
lations such as their implementation in high-order finite difference codes, problems 
with defiltering, and modifications to accomodate integrations to solid boundaries axe 
elaborated on. These as well as other points are discussed in detail along with the 
authors’ views concerning the prospects for future research. 


‘This research was supported by the National Aeronautics and Space Administra- 
tion under NASA Contract No. NASl-18605 while the authors were in residence at 
ICASE, NASA Langley Research Center, Hampton, VA 23665. 
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1. Introductory Remarks 

There is no question that both direct and large-eddy simulations will play an important role 
in turbulence research for many years to come. The position paper of Professor Reynolds 
makes two major points in this regard: (a) direct numerical simulations (DNS) will be 
used primarily to gain a better understanding of basic turbulence physics and to assess the 
performance of models and theories, and (b) large-eddy simulations (LES) will be used as a 
prediction and design tool for complex turbulent flows of scientific and engineering interest. 
Professor Reynolds also alluded to the point that numerical simulations will not supplant 
the need for physical experiments and turbulent closure models but rather will complement 
them. Physical experiments will be necessary to gain a better understanding of the physics 
of turbulent flows whose high Reynolds numbers and/or geometrical complexities make them 
inaccessible to direct numerical simulations. Likewise, turbulent closure models will probably 
need to be used for engineering design purposes in the forseeable future. 

We agree with most of the points made in the Reynolds position paper, but feel that it 
would be useful to clarify and elaborate on some issues concerning the accomplishments and 
limitations of direct and large-eddy simulations. In particular, some clarifications are needed 
in regard to the current status of numerics for turbulence simulations, and Renormalization 
Group Methods. In addition, we shall present our views concerning some critical problems 
that must be overcome before large-eddy simulations can be used on a routine basis for the 
prediction and control of engineering turbulent flows. These points will be discussed in detail 
along with our views of the prospects for future research. 
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2 . Numerics and Mathematical Theory 


The choice of a numerical method for incompressible DNS and LES remains a compromise 
between the accuracy produced by a high-order method and the speed achieved by a low- 
order scheme. By now there is overwhelming numerical evidence that the spatial accuracy 
in direct simulations ought to be at least fourth-order, while there is substantial sentiment 
that large-eddy simulations should be high-order as well. On this last point, however, con- 
siderable philosophical disagreement remains. One camp argues that subgrid-scale modeling 
inaccuracies vitiate any practical advantages of high-order methods over simple second-order 
schemes (Schumann, Grotzbach, and Kleiser [1]). The other camp contends that high-order 
schemes permit clear distinctions to be made between numerical and modeling errors. In 
special cases, such as problems amenable to tensor-product grids, many (including our own), 
but not all, groups have preferred fully spectral schemes (see Canuto, Hussaini, Quarteroni, 
and Zang [2]) to high-order finite-difference or finite-element methods (Browning and Kreiss 
[3]). For more general geometries the choice is between fourth- or perhaps sixth-order finite- 
difference methods (see Krist [4] and Rai and Moin [5]) and spectral domain decomposition 
methods (reviewed in Chapter 4 of Zang, Streett, and Hussaini [6]). The paper by Rai and 
Moin [5], referred to by Reynolds, describes a (mostly) sixth-order finite-difference scheme 
that shows promise, although it has been tested thus far only on a tensor-product grid for 
which fast elliptic solvers are available and for which the deterioration in accuracy near the 
boundary seems not to be significant. We must note, however, that this paper contains no 
evidence to support Reynolds’ speculation that high-order finite-difference methods may be 
more efficient than spectral domain decomposition approaches. It will take several years 
for sufficient experience with both types of methods to accumulate before even tentative 
conclusions about their relative efficiency can be drawn. We suspect that in the end neither 
approach will prove uniformly preferable. 

We concur with Professor Reynolds’ statement that the key to the efficiency of these 
methods is rapid elliptic solvers, and observe that the multigrid community has already pro- 
duced fairly general techniques that are applicable to complex three-dimensional geometries. 
What remains is for these to be implemented in DNS and LES codes (see Krist and Zang 
[7] and Krist [4] for some DNS applications in tensor-product environments). Moreover, 
the prospects are quite high that the computers that achieve the 100 Gigaflops to 1 Teraflop 
speeds projected by Reynolds and many others will be highly parallel machines with program- 
ming considerations quite unlike those for conventional supercomputers (Cray, Cyber-205, 
Fujitsu, etc.) on which most large-scale turbulence simulations have been performed. Thus, 
a major effort is required on highly parallel algorithms for turbulence simulations in general 
geometries. 

Professor Reynolds notes that Horiuti [8] demonstrated that the so-called Arakawa form 
(more properly called the skew-symmetric form) of the nonlinear terms in the momentum 
equation was superior to the rotation form for LES simulations. However, as documented 
for a variety of fully spectral transition and turbulence simulations by Zang [9], the problem 
with the rotation form is not one of greater truncation error but rather one of greater 
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aliasing errors. Indeed, Zang’s results indicate that the difference between aliased and de- 
aliased spectral calculations is very much less when the rotation form is replaced by the 
skew-symmetric form. Moreover, there is now substantial additional numerical evidence 
(Ronquist [10]) that even for spectral domain decomposition methods, the skew-symmetric 
form yields much more accurate results than the rotation form. 

Numerical simulations of turbulence have traditionally used time-discretizations with an 
accuracy ranging between second- and fourth-order, with an implicit treatment of the viscous 
terms and an explicit treatment of the advection terms. The codes are typically operated 
close to the stability limit. This approach is certainly adequate for the lower-order statistics 
but appears to be questionable for the higher-order ones since the smaller scales are treated 
less accurately by the time-discretization. 

The major question facing simulations of compressible flow is what to do about the 
shock waves: resolve them, fit them, or capture them. At sufficiently low Reynolds number 
they can, of course, be resolved. This has been the approach adopted, for example, by 
Passot and Pouquet [11], Lele [12], and Erlebacher, Hussaini, and Kreiss [14]. But there are 
pressing technological issues involving strong shocks whose internal structure is impractical 
to resolve numerically. Some basic issues can be addressed numerically by shock-fitting 
techniques (Zang, Hussaini, and Bushnell [13]), but this approach is not practical for most 
flows. The dilemma facing shock-capturing methods is that the numerical viscosity required 
for capturing the shock may well seriously distort small-scale turbulent features of interest 
near the shock. The model problem studied by Zang, Hussaini, and Bushnell [13] is a good 
test bed for shock-capturing methods and has been used recently to check some modern 
upwind methods (Shu and Osher [16], Zang, Drummond, and Kumar [17]). The recent 
studies of free shear layers by Soetrisno, Eberhard, Riley, and McMurtry [18] are one example 
of the use of shock-capturing in this context. 

In recent years mathematicians have made significant progress in the analysis of the 
relevant range of scales in turbulent flow that has important implications for the numerical 
simulation of turbulence. Foias, Manley, and Temam [19] estimate the cut-off wavenumber 
at which exponential decay of the spectrum ensues. Henshaw, Kreiss, and Reyna [20], [21] 
show that the minimum scale of turbulence is inversely proportional to the square root of the 
Reynolds number based on the viscosity and the maximum of the vorticity. This information 
can be used to predict accurately the resolution requirements for DNS. The latter work also 
shows that in two-dimensional incompressible turbulence, the flow passes from a stage with 
a k ~ 3 spectrum to a stage with a k~ A spectrum. In the first stage the exact dissipation 
mechanism is extremely important. For example, a hyperviscosity subgrid-scale dissipation 
produces a physically incorrect result. In the second stage the enstrophy decays very slowly 
and the Euler equations are a good approximation. Results for three-dimensional turbulence 
appear to require an assumption on the maximum vorticity in the flow. Henshaw, Kreiss, 
and Reyna [20] show rigorously that if the maximum norm of the vorticity scales as i 
then the Kolmogorov power law follows. 
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3. Direct Numerical Simulations 


Incompressible Flows 

Existing data bases on homogeneous turbulence have been quite useful in studying tur- 
bulence structure as discussed in the Reynolds position paper. The most useful data bases 
have been those on plane shear, plane strain, and the axisymmetric contraction/expansion. 
While the homogeneous data bases at NASA Ames/Stanford based on the pioneering work 
of Rogallo [22] are quite extensive, it should be mentioned that there are some notable de- 
ficiencies. For instance, no simulations have been conducted on plane shear or plane strain 
with rotations (only coarsely resolved large-eddy simulations were conducted for rotating 
shear flow by Bardina, Ferziger, and Reynolds [23]). Such flows can be extremely useful in 
calibrating one-point turbulence closure models as demonstrated by Speziale and MacGiolla 
Mhuiris [24]. The large-eddy and direct simulations on rotating isotropic turbulence con- 
ducted by Bardina, Ferziger, and Rogallo [25] and Speziale, Mansour, and Rogallo [26] are 
interesting but possibly inappropriate for calibrating one-point turbulence closure models; 
the suppression of the energy cascade by a system rotation (which is the primary physi- 
cal effect manifested in rotating isotropic turbulence), is decidedly a two-point phenomenon 
that is not likely to be described properly within the framework of the usual Reynolds stress 
closures. 

The DNS of inhomogeneous turbulent flows have primarily concentrated on boundary 
layer and plane and curved channel flow. As pointed out in the Reynolds position paper, 
this work has shed some new light on the turbulence structure near solid boundaries. How- 
ever, a word of caution is but proper. The Reynolds number of the channel flow simulations 
are so low that serious questions can be raised about the applicability of the results to high 
Reynolds number turbulent channel flow with fully developed turbulence statistics. A mini- 
mization of Reynolds number effects (i.e., for the turbulence statistics outside the wall layer 
to be independent of the Reynolds number) would require at least an eight-fold increase in 
Reynolds number (i.e., from Re ss 4000 to Re « 32,000). Such computations would require 
approximately a three orders of magnitude increase in computational power - a figure that 
will be achievable in the not too distant future. Such an increase in computational capac- 
ity would also allow for homogeneous turbulence simulations that are within the Reynolds 
number range of many of the physical experiments that have been conducted. Hence, while 
we feel that some significant strides have been made in studying turbulence physics based on 
DNS data bases, we believe that results of a lasting impact are more likely to be achieved by 
the higher Reynolds number data bases that are likely to be generated after the next decade 
(i.e., by the turn of the century). 

The work of Gilbert and Kleiser [27] has shown that the constant flow rate formulation of 
turbulent channel flow is more economical than the constant pressure gradient formulation 
that has been customary. The recent landmark calculation by Gilbert [28], in collaboration 
with Kleiser, of the complete transition to turbulence process in channel flow bears mention. 
For the first time they have reliably computed both the transitional region and the fully 
turbulent regime. As recent calculations by Zang and Krist [29] suggest, this achievement 
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was due in part to a fortunate choice of parameters, for at higher Reynolds numbers or 
with different initial conditions the transitional regime may actually require substantially 
more resolution than the turbulent regime. Nevertheless, Gilbert’s data does offer the first 
opportunity to examine the final stages of breakdown into turbulence. We should also note 
that the results of Zang and Krist for transitional flow clearly indicate the mechanisms 
responsible for inverted vortices. 

Compressible Flows 

Work done at Langley Research Center (Erlebacher, Hussiani, and Kreiss [14] and Speziale, 
Erlebacher, Zang, and Hussaini [15]) has also shown that in compressible DNS of homoge- 
neous turbulent flow, the flow statistics and the kinetic energy remain mostly incompressible 
when starting from essentially incompressible initial conditions. We have also shown that this 
is primarily due to the initial conditions for 2-D isotropic turbulence. Passot and Pouquet 
[11] state that it is the relative importance of the fluctuating Mach number, Ma , relative 
to the rms density fluctuations, p rma which determines the final outcome of the turbulent 
structures, and these results have been confirmed by Erlebacher, et al.[14]. If Ma 2 « 1 and 
Prma = 0(1), then weak shocks prevail and most of the kinetic energy is compressible; on 
the other hand, if Ma 2 << 1 and p rma — O(Ma), no shocks occur and most of the kinetic 
energy is solenoidal. The energy spectra and contours of the divergence of velocity associated 
with these two regimes which are taken from the direct simulations of Erlebacher et al. [14] 
on 2-D isotropic turbulence are illustrated in Figures 1-4. Of course, strong shocks appear 
if Ma = 0(1). The curvature of these shocks then generates localized regions of intense 
vorticity. Passot and Pouquet as well as Erlebacher, Hussaini and Kreiss [14] have observed 
shocks/shocklets in 3-D isotropic, homogeneous turbulence as well. The dynamics seem to 
agree with the 2-D case, but at a higher Ma in 3-D. 

The fact that shocks have not yet been found in three-dimensional shear flow simula- 
tions at NASA Ames Research Center does not prove that they do not exist. Erlebacher, 
et al have already found eddy shocklets in three-dimensional isotropic turbulence, and it 
seems reasonable to expect that a proper choice of initial conditions will induce them in 
shear layers as well. There is a large parameter regime to study. Although not all of the 
important mechanisms of supersonic compressible turbulence can be extracted from two- 
dimensional simulations, the uniquely compressible effects arising from the interactions of 
large-scale shocks with turbulence are present in two-dimensions and can provide some useful 
information. 


4. Large-Eddy Simulations 

It has been speculated since the mid 1970’s that large-eddy simulations would some day 
become a routine tool for the analysis of complex turbulent flows that cannot be predicted 
by existing turbulence models. To date, the accomplishments of LES have been somewhat 
disappointing. As mentioned in the position paper of Reynolds, LES has “paved the way 
for the use of direct numerical simulations in turbulence”. However, there have been few 


5 


practical uses of LES outside of the geophysical fluid dynamics community (e.g., in weather 
prediction) and there the details of the turbulence structure do not seem to be so important. 
There is no doubt that the lack of wide availability of supercomputers and the redirection of 
research efforts toward direct simulations have stunted the growth of LES as pointed out in 
Reynolds paper. However, it is our opinion that there are some major difficulties with LES 
that need to be overcome before it can yield reliable and economically feasible predictions 
for the complex turbulent flows of scientific and engineering interest. These problems are as 
follows: 

(i) the implementation of LES in spectral domain decomposition or high-order finite- 
difference codes so that complex geometries can be treated; 

(ii) the development of improved subgrid scale models for strongly inhomogeneous 
turbulent flows (e.g., flows with localized regions of relaminarization or large mean 
velocity gradients); 

(iii) the development of reliable a priori tests for the screening of new subgrid scale 
models; 

(iv) the problem of defiltering; 

(v) the problem of modifying subgrid scale models to accomodate integrations to solid 
boundaries. 

In so far as point (i) is concerned, it must be emphasized that most of the successes of 
LES have been for either homogeneous turbulent flows or simple inhomogeneous turbulent 
flows (e.g., channel flow). To accommodate complex geometries, either spectral domain 
decomposition techniques or finite- difference or finite-element methods must be used. Proper 
filtering techniques which are appropriate for such methods must be developed. Furthermore, 
to compute turbulence statistics, either ensemble or time averages must be implemented 
in lieu of the spatial averages used for homogeneous turbulent flows or the simple planar 
averages used in turbulent channel flow. These complicating features can increase by up to 
two orders of magnitude the computational time needed in comparison to the LES of more 
simple homogeneous turbulent flows. In fact, the latter problem of constructing ensemble or 
time averages instead of spatial averages makes it, in our opinion, highly unlikely that the 
LES of complex three-dimensional flows will be less expensive than computations with more 
sophisticated Reynolds stress models such as second-order closures. 

In order to obtain reliable predictions for strongly inhomogeneous turbulent flows with 
large mean field gradients there is the likely prospect that simple eddy viscosity models will 
not suffice as alluded to in point (ii). Subgrid scale transport models that account for non- 
local and history-dependent effects may well be needed (cf., Deardorff [30]). Such models 
have the potential of increasing the computational efforts required by as much as a factor 
of five. Furthermore, there is the need for a substantial research effort to develop transport 
models that can yield reliable predictions and are well-behaved. 
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The more efficient development of improved subgrid scale models would be aided con- 
siderably by the development of reliable a priori tests as mentioned in point (iii). The 
commonly used a priori tests (see Clark, Ferziger, and Reynolds [31] and Bardina, Ferziger, 
and Reynolds [23]) are rather unreliable in predicting the performance of LES as alluded to 
in the position paper of Reynolds. To elaborate further on this point, it has long been spec- 
ulated that the success of LES hinges on the use of a subgrid scale model that dissipates the 
right amount of energy in order to compensate for the energy cascade to scales that are not 
resolved by the computational mesh. However, the Smagorinsky model has been shown to 
yield a correlation coefficient of roughly 0.5 based on dissipation. To gain an appreciation as 
to how poor this correlation is, it should be noted that the correlation between the functions 
y = x and y = exp(— x) ( two functions with entirely different qualitative behavior) is 0.696 
on the interval [0,1]! Yet, despite this extremely poor correlation of 0.5, the Smagorinsky 
model has performed well in many large-eddy simulations. Consequently, it is clear that 
alternative means of correlating subgrid scale models with the results of DNS are needed. 
One possible alternative might be to base the correlations on dissipation subject to the con- 
dition that the subgrid scale dissipation is some specified, substantial amount larger than 
the viscous dissipation. This issue needs to be explored in more detail in the near future. 


Anisotropy Tensor 

Large-Eddy Simulations 

Experiments 

bn 

0.305 

0.201 

b 2 2 

-0.265 

-0.148 

b\2 

-0.15 

-0.142 


Table 1. Comparison of the anisotropy tensor obtained from the large-eddy simula- 
tion of Bardina, et al. [23] for homogeneous turbulent shear flow with the physical 
experiments of Tavoularis and Corrsin [32]. 


If higher-order turbulence statistics are needed beyond the mean velocity, the problem of 
defiltering arises (point (iv)). This problem becomes most critical when more than 10-20% 
of the turbulent kinetic energy is in the subgrid scale motion. Typically as much as 50% 
of the turbulent kinetic energy is in the SGS fields in practical LES. (Certainly if the value 
were only 10%, little computational savings would be gained for a given Reynolds number!). 
With so much energy in the subgrid scale motions, considerable errors can be introduced in 
the computation of second and higher-order moments (see Table 1). This points to the need 
for defiltering wherein an estimate must be made for the contributions from the subgrid scale 
fields (see Bardina et al. [23]). It must be said at the outset that there is no general, reliable 
method for defiltering. The problem of defiltering is tantamount to solving a Fredholm 
integral equation of the first kind 

¥=y G(x-x / )$(x , )d 3 x' ) (1) 
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where <f> is any given flow variable, $ is its filtered form, and G is the filter function. Since a 
Fredholm integral equation is ill-posed in the sense of Hadamard (i.e., any numerical solution 
for $ yields complete uncertainty in regard to the high frequency Fourier components of 4>), 
the higher-order statistics associated with $ can only be crudely estimated. Hence, LES can 
only generally give accurate information on first-order moments. More formal probabilistic 
methods for estimating the SGS contributions to second-order moments should be studied 
in the future. 

Finally, we would like to make a few remarks on point (v) concerning the use of subgrid 
scale models near a solid boundary. Most of the LES conducted have used ad hoc modifica- 
tions such as Van Driest damping near solid boundaries (e.g., Schumann [33], and Moin and 
Kim [34]). While such an approach is satisfactory for attached flows in simple geometries, it 
cannot be reliably applied to turbulent flows with separation or complex geometrical singu- 
larities. Here, we believe the Renormalization Group (RNG) approach of Yakhot and Orszag 
[35] holds more promise than other existing methods in bridging the eddy viscosity to the 
molecular viscosity as a solid boundary is approached. The LES of turbulent channel flow 
based on RNG (Orszag 1987, private communication) shows considerable promise in pre- 
dicting the turbulence structure near solid boundaries. In our opinion, the RNG approach 
to LES deserves more future research. 


5. Concluding Remarks 

As outlined in the Reynolds position paper, significant strides have been made in the analysis 
of turbulence physics based on direct numerical simulations of the Navier-Stokes equations. 
The low Reynolds numbers of these simulations do somewhat minimize the long-range im- 
pact that the specific results obtained will have on our understanding of basic turbulence 
physics. Nevertheless, an important methodology has been developed which, by the turn of 
the century (with a projected three orders of magnitude increase in computational speed), 
will allow for the direct simulation of basic turbulent flows at Reynolds numbers that are 
high enough to be representative of the physical mechanisms observed in the turbulent flows 
of scientific and engineering interest. 

In so far as large-eddy simulations are concerned, we believe that there are major op- 
erational problems that must be overcome (as outlined in this discussion) before they can 
become a more commonly used tool for the analysis of complex engineering flows. Nonethe- 
less, we feel that there is the potential for making a major impact with LES in scientific and 
engineering computations so that research along these lines should be vigorously pursued. 
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Figure 1. Energy spectra for 2-D isotropic decay of compressible 
turbulence taken from the direct simulations of Erlebacher, et 
al. [14]: the eddy shocklet regime; Re A = 10, Ma = 0.028, 
Prmi =0.10. 
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Figure 2. Contours of the divergence of velocity for 2-D isotropic decay 
of compressible turbulence taken from the direct simulations 
of Erlebacher, et al. [14]: the eddy shocklet regime; Re* = 10, 
Ma = 0.028, Prmt = 0.10. 
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Figure 3. Energy spectra for 2-D isotropic decay of compressible 
turbulence taken from the direct simulations of Erlebacher, et 
al. [14]: the quasi-incompressible regime;. Re A = 10, Ma = 0.028, 
Prmt = 0. 
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Figure 4. 
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Contours of the divergence of velocity for 2-D isotropic decay 
of compressible turbulence taken from the direct simulations 
of Erlebacher, et al. [14]: the quasi-incompressible regime; 


Re A = 10, Ma = 0.028, p rm . = 0. 
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